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Abstract 

Many of the static and dynamic properties of an atomic Bose-Einstein condensate (BEC) are usually 
studied by solving the mean-field Gross-Pitaevskii (GP) equation, which is a nonlinear partial differential 
equation for short-range atomic interaction. More recently, BEG of atoms with long-range dipolar atomic 
interaction are used in theoretical and experimental studies. For dipolar atomic interaction, the GP equation 
is a partial integro-differential equation, requiring complex algorithm for its numerical solution. Here we 
present numerical algorithms for both stationary and non-stationary solutions of the full three-dimensional 
(3D) GP equation for a dipolar BEG, including the contact interaction. We also consider the simplified 
one- (ID) and two-dimensional (2D) GP equations satisfied by cigar- and disk-shaped dipolar BECs. We 
employ the split-step Crank-Nicolson method with real- and imaginary-time propagations, respectively, for 
the numerical solution of the GP equation for dynamic and static properties of a dipolar BEC. The atoms 
are considered to be polarized along the z axis and we consider ten different cases, e.g., stationary and non- 
stationary solutions of the GP equation for a dipolar BEC in ID (along x and z axes), 2D (in x — y and x — z 
planes), and 3D, and we provide working codes in Fortran 90/95 and C for these ten cases (twenty programs 
in all). We present numerical results for energy, chemical potential, root-mean-square sizes and density of 
the dipolar BECs and, where available, compare them with results of other authors and of variational and 
Thomas-Fermi approximations. 
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Program summary 

Program title: (i) imagldZ, (ii) imagldX, (iii) iinag2dXY, (iv) imag2dXZ, (v) imagSd, (vi) realldZ, (vii) 
realldX, (viii) real2dXY, (ix) real2dXZ, (x) real3d 
Catalogue identifier. AEWL_vl_0 

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEWL_vl_0.htinl 
Program obtainable from: CPC Program Library, Queens University, Belfast, N. Ireland 
Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/hcence/licence.html 
Distribution format: tar.gz 
Programming language: Eortran 90/95 and C 

Computer: Any modern computer with Fortran 90/95 or C language compiler installed. 

Operating system: Linux, Unix 

RAM: 1 GB (i,ii), 2 GB (iii,iv), 4 GB (v), 2 GB (vi,vii), 4 GB (viii,ix), 8 GB (x) 

Classification: 2.9, 4.3, 4.12 

Nature of problem: These programs are designed to solve the time-dependent nonlinear partial differential 
Gross-Pitaevskii (GP) equation with contact and dipolar interactions in one, two or three space dimensions 
in a harmonic anisotropic trap. The GP equation describes the properties of a dilute trapped Bose-Einstein 
condensate. 

Solution method: The time-dependent GP equation is solved by the split-step Grank-Nicolson method by 
discretizing in space and time. The discretized equation is then solved by propagation, in either imaginary 
or real time, over small time steps. The contribution of the dipolar interaction is evaluated by a Fourier 
transformation (FT) to momentum space using a convolution theorem. The method yields the solution of 
stationary and/or non-stationary problems. 

Additional comments: This package consists of 12 programs, see Program title, above. Fortran 90/95 and 
G versions are provided for each of the 6 programs. For the particular purpose of each program please see 
below. 

Running time: Minutes on a medium PC (i, ii), tens of minutes on a medium PC (iii, iv), tens of minutes 
on a good workstation (v, vi). 
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1. Introduction 


After the experimental realization of atomic Bose-Einstein condensate (BEC) of alkali-metal and some 
other atoms, there has been a great deal of theoretical activity in studying the statics and dynamics of 
the condensate using the mean-field time-dependent Gross-Pitaevskii (GP) equation under different trap 
symmetries [1]. The GP equation in three dimensions (3D) is a nonlinear partial differential equation in 
three space variables and a time variable and its numerical solution is indeed a difficult task specially for 
large nonlinearities encountered in realistic experimental situations [2]. Very special numerical algorithms 
are necessary for its precise numerical solution. In the case of alkali-metal atoms the atomic interaction 
in dilute BEG is essentially of short-range in nature and is approximated by a contact interaction and at 
zero tempearture is parametrized by a single parameter in a dilute BEC — the s-wave atomic scattering 
length. Under this approximation the atomic interaction is represented by a cubic nonlinearity in the GP 
equation. Recently, we published the Fortran [3] and C [3] versions of useful programs for the numerical 
solution of the time-dependent GP equation with cubic nonlinearity under different trap symmetries using 
split-step Crank-Nicolson scheme and real- and imaginary-time propagations. Since then, these programs 
enjoyed widespread use [5]. 

More recently, there has been experimental observation of BEC of ^^Cr i, 164Dy [7] and i®®Er i atoms 
with large magnetic dipole moments. In this paper, for all trap symmetries the dipolar atoms are considered 
to be polarized along the z axis. In these cases the atomic interaction has a long-range dipolar counterpart in 
addition to the usual contact interaction. The s-wave contact interaction is local and spherically symmetric, 
whereas the dipolar interaction acting in all partial waves is nonlocal and asymmetric. The resulting GP 
equation in this case is a partial integro-differential equation and special algorithms are required for its 
numerical solution. Different approaches to the numerical solution of the dipolar GP equation have been 
suggested [iiiniiniiiiiiiiiii. Yi and You m solve the dipolar GP equation for axially-symmetric trap 
while they perform the angular integral of the dipolar term, thus reducing it to one in axial (z,z') and 
radial (p, p') variables involving standard Elliptical integrals. The dipolar term is regularized by a cut-off at 
small distances and then evaluated numerically. The dipolar GP equation is then solved by imaginary-time 
propagation. Goral and Santos El treat the dipolar term by a convolution theorem without approximation, 
thus transforming it to an inverse Fourier transformation (FT) of a product of the FT of the dipolar potential 
and the condensate density. The FT and inverse FT are then numerically evaluated by standard fast Fourier 
transformation (FFT) routines in Cartesian coordinates. The ground state of the system is obtained by 
employing a standard split-operator technique in imaginary time. This approach is used by some others 
m- Ronen et al. perform the angular integral in the dipolar term using axial symmetry. To evaluate 
it, in stead of FT in x,y, and z El, they use Hankel transformation in the radial p variable and FT in 
the axial z variable. The ground state wave function is then obtained by imaginary-time propagation and 
dynamics by real-time propagation. This approach is also used by some others m- Bao et al. use Euler sine 
pseudospectral method for computing the ground states and a time-splitting sine pseudospectral method for 
computing the dynamics of dipolar BECs [9]. Blakie et al. solve the projected dipolar GP equation using 
a Hermite polynomial-based spectral representation m- Lahaye et al. use FT in x,y, and z to evaluate 
the dipolar term and employ imaginary- and real-time propagation after Crank-Nicolson discretization for 
stationary and nonstationary solution of the dipolar GP equation El- 

Here here we provide Fortran and C versions of programs for the solution of the dipolar GP equation in a 
fully anisotropic 3D trap by real- and imaginary-time propagation. We use split-step Crank-Nicolson scheme 
for the nondipolar part as in Refs. mm and the dipolar term is treated by FT in x, y, z variables. We also 
present the Fortran and C programs for reduced dipolar GP equation in one (ID) and two dimensions (2D) 
appropriate for a cigar- and disk-shaped BEC under tight radial (p) and axial (z) trapping, respectively 
El- In the ID case, we consider two possibilities: the ID BEC could be aligned along the polarization 
direction z or be aligned perpendicular to the polarization direction along x axis. Similarly, in the 2D case, 
two possibilities are considered taking the 2D plane as x — y, perpendicular to polarization direction z or 
as a: — z containing the polarization direction. This amounts to five different trapping possibilities — two in 
ID and 2D each and one in 3D — and two solution schemes involving real- and imaginary-time propagation 
resulting in ten programs each in Fortran and C. 
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In Sec. [^we present the 3D dipolar GP equation in an anisotropic trap. In addition to presenting the 
mean-field model and a general scheme for its numerical solution in Secs. 2.1 and 2.2 we also present two 
approximate solution schemes in Secs. |2.3| and |2.4[ — Gaussian variational approximation and Thomas- 
Fermi (TF) approximation — in this case. The reduced ID and 2D GP equations appropriate for a cigar- 
and a disk-shaped dipolar BEG are next presented in Secs. |2.5| and |2.6[ respectively. The details about 
the computer programs, and their input/output files, etc. are given in Sec. The numerical method and 
results are given in Sec. Finally, a brief summary is given in Sec. 


2. Gross-Pitaevskii (GP) equation for dipolar condensates in three dimensions 

2.1. The mean-field Gross-Pitaevskii equation 

At ultra-low temperatures the properties of a dipolar condensate of iVat atoms, each of mass m, can be 
described by the mean-field GP equation with nonlocal nonlinearity of the form: nnum 


ih 


dffivG) 

dt 


v72 , T/ ^ ^ 


|.^(r,i)|^ + Wat f (/dd(r - r') |0(r',f)t dr' 


( 1 ) 


where J dr\(j){r,t)\'^ = 1. The trapping potential, Vtrap is assumed to be fully asymmetric of the form 

Ptrap(r) = + ujIz^) 

where uJx, 0 Jy and ujz are the trap frequencies, a the atomic scattering length. The dipolar interaction, for 
magnetic dipoles, is given by [TTJ ITH] 


C/dd(R) 


1 — 3 cos^ 8 
47r |RP 


( 2 ) 


where R = r — r' determines the relative position of dipoles and 9 is the angle between R and the direction 
of polarization z, /rg is the permeability of free space and p, is the dipole moment of the condensate atom. 
To compare the contact and dipolar interactions, often it is useful to introduce the length scale Odd = 
fj.op^m/{l2TTh'^) [5]. 

Convenient dimensionless parameters can be defined in terms of a reference frequency ui and the corre¬ 
sponding oscillator length I = Using dimensionless variables r = r/Z,R = R/Z,a = a/l,ddd = 

= tu!, X = x/l,y = y/l,z = z/l.f) = Eq. Q can be rewritten (after removing the overhead 

bar from all the variables) as 


.d(t){r,t) 
* dt 


+ 47raA^at|<^P + SlVatOdd 


Vif{RM{r\t)\^dP 


(^(r,t), (3) 


with 


KiT(R) 


1 — 3 cos^ 9 

|RP 


( 4 ) 


where 7 = u!xfdj,v = u}y/uj,\ = Wz/w. The reference frequency oj can be taken as one of the frequencies 
(jjx,ojy or ujz or their geometric mean (uJxUJyUJz)^'^^■ In the following we shall use Eq. ^ where we have 
removed the ‘bar’ from all variables. 

Although we are mostly interested in the numerical solution of Eq. ([^, in the following we describe 
two analytical approximation methods for its solution in the axially-symmetric case. These approximation 
methods — the Gaussian variational and TF approximations — provide reasonably accurate results under 
some limiting conditions and will be used for comparison with the numerical results. Also, we present 
reduced ID and 2D mean-field GP equations appropriate for the description of a cigar and disk-shaped 
dipolar BEG under appropriate trapping condition. The numerical solution and variational approximation 
of these reduced equations will be discussed in this paper. A brief algebraic description of these topics are 
presented for the sake of completeness as appropriate for this study. For a full description of the same the 
reader is referred to the original publications. 
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2.2. Methodology 

We perform numerical simulation of the 3D GP equation (|^ using the split-step Crank-Nicolson method 
described in detail in Ref. [3]. Here we present the procedure to include the dipolar term in that algorithm. 
The inclusion of the dipolar integral term in the GP equation in coordinate space is not straightforward due 
to the singular behavior of the dipolar potential at short distances. It is interesting to note that this integral 
is well defined and finite. This problem has been tackled by evaluating the dipolar term in the momentum 
(k) space, where we do not face a singular behavior. The integral can be simplified in Fourier space by 
means of convolution as 

I dr'Vl^ir - v')n{v\t) = J t), (5) 

where n(r,t) = |^(r,t)p. The Fourier transformation (FT) and inverse FT, respectively, are defined by 


H(k) = I drH(r)e»'^-, A{r) = ^ J dkH(k)e-*'^ G (6) 

The FT of the dipole potential can be obtained analytically [13] 

i5i>(k)sf(.3D(k)=|(S_l), (7, 

SO that 

J dr'Vif{r-r')n{r',t) = ^ J ^^^e-*'" '’h3D(k)n(k, t). (8) 

To obtain Eq. Q, first the angular integration is performed. Then a cut-off at small r is introduced to 
perform the radial integration and eventually the zero cut-off limit is taken in the final result as shown 
in Appendix A of Ref. [19]. The FT of density n(r) is evaluated numerically by means of a standard 
FFT algorithm. The dipolar integral in Eq. ([^ involving the FT of density multiplied by FT of dipolar 
interaction is evaluated by the convolution theorem ([^ . The inverse FT is taken by means of the standard 
FFT algorithm. The FFT algorithm is carried out in Gartesian coordinates and hence the GP equation is 
solved in 3D irrespective of the symmetry of the trapping potential. The dipolar interaction integrals in ID 
and 2D GP equations are also evaluated in momentum spaces. The solution algorithm of the GP equation 
by the split-step Grank-Nicolson method is adopted from Refs. Ell]. 

The 3D GP equation is numerically the most difficult to solve involving large RAM and GPU time. 
A requirement for the success of the split-step Grank-Nicolson method using a FT continuous at the origin 
is that on the boundary of the space discretization region the wave function and the interaction term should 
vanish. For the long-range dipolar potential this is not true and the FT Q is discontinuous at the origin. 
The space domain (from —oo to -l-oo) cannot be restricted to a small region in space just covering the spatial 
extention of the BEG as the same domain is also used to calculate the FT and inverse FT used in treating 
the long-range dipolar potential. The use and success of FFT implies a set of noninteracting 3D periodic 
lattice of BEGs in different unit cells. This is not true for long-range dipolar interaction which will lead 
to an interaction between BEGs in different cells. Thus, boundary effects can play a role when finding the 
FT. Hence a suffciently large space domain is to be used to have accurate values of the FT involving the 
long-range dipolar potential. It was suggested [12] that this could be avoided by truncating the dipolar 
interaction conveniently at large distances r = R so that it does not affect the boundary, provided R is 
taken to be larger than the size of the condensate. Then the truncated dipolar potential will cover the whole 
condensate wave function and will have a continuous FT at the origin. This will improve the accuracy of a 
calculation using a small space domain. The FT of the dipolar potential truncated at r = i?, as suggested 
in Ref. [12] . is used in the numerical routines 


Kir(k) 


f(f-0 


^cos(fci?) ^sin(fcR) 
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A: = |k|. 


(9) 









Needless to say, the difFculty in using a large space domain is the most severe in 3D. In 3D programs the 
cut-off R of Eq. (9) improves the accuracy of calculation and a smaller space region can be used in numerical 
treatment. In ID and 2D, a larger space domain can be used relatively easily and no cut-off has been used. 
Also, no convenient and effcient analytic cut-off is known in ID and 2D [T^. The truncated dipolar potential 
0 has only been used in the numerical programs in 3D, e.g., imag3d* and real3d*. In all other numerical 
programs in ID and 2D, and in all analytic results reported in the following the untruncated potential Q 
has been used. 


2.3. Gaussian variational approximation 

In the axially-symmetric case (j = i/), convenient analytic Lagrangian variational approximation of Eq. 
([^ can be obtained with the following Gaussian ansatz for the wave function |20] 


(/>(r,t) = 


,-3/4 


Wp{t)yJW:,{t) L 2 ^ 2 (t) 2wl{t) 


+ ia{t)p'^ + if3{t)z''^ 


( 10 ) 


where r = {p, z}, p = {a;, y}, Wp{t) and Wz{t) are widths and all:.) and /3(t) are chirps. The time dependence 
of the variational parameters Wp{t), Wz{t), a{t) and (3(t) will not be explicitly shown in the following. 

The Lagrangian density corresponding to Eq. 0 is given by 


/:(r) = 2 


(l)(r) 


dt 


- ^*{^) 


d(t){-r) 

dt 


+ ^ + 27raiVat|<(>(r)|‘‘ 


3Udd A^at 


|</.(r)p / Edd(R)|<(>(r')Pdr'. 


( 11 ) 


Consequently, the effective Lagrangian L = J £(r)dr (per particle) becomes [51ET] 

1 


r 2 . , 7 Wp 

L-u^a+ — + — 


11 I 9 2 2 I 2o2| -Aat[« add/(K)] 

- ^ + 2.U} a +i^zP d 

4 2ujI 4a;2 ^ 

H ^ 


( 12 ) 


The Euler-Lagrangian equations with this Lagrangian leads to the following set of coupled ordinary differ¬ 
ential equations (ODE) for the widths Wp and Wz [22] 

A^at [ 2 a - addff(K)] 


2 1 

Wp+-1 Wp = ^ + 


Wz + X Wz = ^ + 


1 2A^at [a - addc(K)] 


W^Wz 


wi 




w^wl 


with K = Wpjwz and 


, , 2 — 7k^ — 4k^ -|- 9K^d{K) 

=-(ITTpP-. 


z{k) = 


1 + IOk^ — 2 k^ — 9K^d{K) 

( 1 -k 2)2 ^ 


... 1 - 1 - 2^2 — Sn^dlK) ,, , atanhVl — ^2 

/(k) = , d{K) = 


1 — K? 




(13) 

(14) 

(15) 

(16) 
(17) 


The widths of a (time-independent) stationary state are obtained from Eqs. and ( [T4| ) by setting 

Wp = Wz = 0. The energy (per particle) of the stationary state is the Lagrangian (12) with a = /3 = 0, e.g.. 


E 


2?«2 4x(;2 


Natja - addfin)] 1 wl ^ X^wl 

\/2ttWzw‘1 2 4 


The chemical potential p = dE/dNg,t of the stationary state is given by [22] 

11 2Afat[a-add/(«:)] , 


^ 2w'^ 4r(;2 


'J2/K' 


TTWzW^ 


(18) 


(19) 
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2.4- Thomas-Fermi (TF) approximation 

In the time-dependent axially-symmetric GP equation ([^, when the atomic interaction term is large 
compared to the kinetic energy gradient term, the kinetic energy can be neglected and the useful TF 
approximation emerges. We assume the normalized density of the dipolar BEG of the form [H |231 [2H [21] 


n(r,t) = 


8TrRl{t)R,{t) 




( 20 ) 


where Rp{t) and Rz(t) are the radial and axial sizes. The time dependence of these sizes will not be explicitly 
shown in the following. Using the parabolic density (20), the energy functional Etf may be written as [24] 


Etf = F\ 


trap 


E\nt — 


N{2'y^RI + X^Rl) 

14 


15 47raiV^.t 
^ R^R^ 




( 21 ) 


where k = Rp/Rz is the ratio of the condensate sizes and /(k) is given by Eq. ( [Tt] ). In Eq. (21), Etrap is 
the energy in the trap and Uint is the interaction or release energy in the TF approximation. In the TF 
regime one has the following set of coupled ODEs for the evolution of the condensate sizes [2H] : 


Rp = -RpT 

Rz = -X'^Rz 


I5aiVat 

RpRz 

15a7Vat 


i ?2 


Rj 


Rl 


F 


^dd 

a 

2 add 

a 


Rl 

I 

m 


3 /(k) 
2Rj-Rl,, 
, 3 /(k) 


2RI-RI 


( 22 ) 

(23) 


The sizes of a stationary state can be calculated from Eqs. (22) and (23) by setting the time derivatives Rp 
and Rz to zero leading to the transcendental equation for k |23j 


3/v 


2 ^dd 


i-k 


2r 


/(^) 

1 — k^ 


- I 


+ (^-1 

V a 


k^- 


A' 


= 0 , 


(24) 


and 


Rp = 


l5aNa.tK 

rf 


1 -f 


Odd ( 3 K^/(/c) 
2 1 - k2 


- 1 


-I 1/5 


(25) 


with Rz = Rp/K. The chemical potential is given by 


h'TF — .^trap T — 


15 47roA^a 
Stt R/Rz 


, Odd rf-\ 

1 -/('«) 

a 


(26) 


We have the identities ExF/Nat = 5/tT^’/7, Ei„t/iVat = 2fiTF/7, Et^^p/N^t = 2 >p.tf/T- 


2.5. One-dimensional GP equation for a cigar-shaped dipolar BEG 
2.5.1. z direction 

For a cigar-shaped dipolar BEG with a strong axially-symmetric (y = 7 ) radial trap (A < iz, 7 ), we assume 
that the dynamics of the BEG in the radial direction is confined in the radial ground state I11IIM1I27] 

(j){p) = exp(-pV2dp)/(dpV^), 7dp = 1, P= {x,y), (27) 

of the transverse trap and the wave function 4>(r) can be written as 


(j){r,t) = (/loizR) X (j){p) 




(t>iD{z,t), 


(28) 
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where 4>iD[z,t) is the effective ID wave function for the axial dynamics and dp is the radial harmonic 
oscillator length. 

To derive the effective ID equation for the cigar-shaped dipolar BEC, we substitute the ansatz ( p^ in 
Eq. multiply by the ground-state wave function (t>{p) and integrate over p to get the ID equation 


* dt 


dl 2afVat|(/>iD| 




Kid (^) = 


27r 

\/2dp 


2 2 

+ 2 i/wi - V^(l + 2w)e^{l - erf(i/w)} 

O 


/ oo 

Kid^{K - z'\)\(j)iD{z',t)\'^dz’ 

-OO 


2.5 


and 


2.6 


(j)iD{z,t), (29) 
(30) 

we use the untruncated 


where w = [Z/{'j2dp)Y, Z = |z —z'|. Here and in all reductions in Secs, 
dipolar potential Q and not the truncated potential The integral term in the ID GP equation ( [2^ is 
conveniently evaluated in momentum space using the following convolution identity 




27r 


kzdp 


where 


/ OO 

-OO 


hioiC) = 


1 


( 27 r) 


dko 


k 2 


- 1 


Rkp)P = 


1 


2'Kd'i 


^^2 1 1,2 

“T 


pOO 

/ du 

Jo 




C = 


kzdp 


(31) 

(32) 

(33) 

(34) 


the following Gaussian ansatz for the wave function 


The ID GP equation p9| can be solved analytically using the Lagrangian variational formalism with 

(35) 


7r-i/4 

(t>iD{z,t) = ! ^ ^ exp 




2wl{t) 




where Wzit) is the width and j3{t) is the chirp. The Lagrangian variational formalism leads to the following 
equation for the width Wz{t) 


Wz{t) + X^Wz[t) = 


1 


+ 


2iVat [a - addc(K)] 


wl{t) 


Wz{t) ' 


(36) 


The time-independent width of a stationary state can be obtained from Eq. (36) by setting Wz{t) = 0. The 
variational chemical potential for the stationary state is given by [ 22 ] 


_ 2iVat[a - addfjk)] 

^ '/^Wzdl 4 


The energy per particle is given by 


E _ Natja - addfjk)] 

Nut 4^2 \/^Wzdj, 4 


(37) 


(38) 


2.5.2. X direction 

For a cigar-shaped dipolar BEC with a strong axially-symmetric {v = A) radial trap (7 < v, X), we assume 
that the dynamics of the BEC in the radial direction is confined in the radial ground state |22l [23 [27] 


cj){p) = exp{-p'^/2dl)/{dpy/X), vd? = 1 , p= {y, z), 


( 39 ) 





































of the transverse trap and the wave function 4>(r) can be written as 






(t>lD{x, t), 


(40) 


where (jjijy^x, t) is the effective ID wave function for the dynamics along x axis and dp is the radial harmonic 
oscillator length. 

To derive the effective ID equation for the cigar-shaped dipolar BEC, we substitute the ansatz (40) in 
Eq. ([^, multiply by the ground-state wave function ^(p) and integrate over p to get the ID equation 


.d(j)iD{x,t) 

* dt 


ilD{Tx) = 


dl 2aNg.t\4>iD\ 


" L 2 ^ 

= dpkxj 72 and 

2 ^ 

dl 

1 f 

ikf 

l^(kp)P 

- (2,r j 

k 2 


/ °° dk 

-or, XTT 


(41) 


72 

2'Kdn 


dxyC '^»h2D(r), 


dpky 

"^=77^ 


T = ^/T^+T^., 


^ 2 d(t) = [2 - 37^6^V{1 - erf(T)}]. 


(42) 

(43) 


To derive Eq. ( |4l] ), the dipolar term in Eq. @ is first written in momentum space using Eq. ([^ and the 
integrations over ky and k^ are performed in the dipolar term. 


2.6. Two-dimensional GP equation for a disk-shaped dipolar BEC 
2.6.1. X — y plane 

For an axially-symmetric (^ = 7 ) disk-shaped dipolar BEC with a strong axial trap (A > v, 7 ), we 
assume that the dynamics of the BEC in the axial direction is confined in the axial ground state 


(j){z) = exp(-zV 2 d^)/( 7 rd 2 )i/ 4 ^ ^ 7VW> 


and we have for the wave function 


(j){r) = (j){z) X (t) 2 D{p,t) 


(7rd2)i/4 



4>2D{P,t), 


(44) 


(45) 


where p = {x,y), (/> 2 d(p, t) is the effective 2D wave function for the radial dynamics and dz is the axial 
harmonic oscillator length. To derive the effective 2D equation for the disk-shaped dipolar BEC, we use 
ansatz ( |4^ i n Eq. ([^, multiply by the ground-state wave function (f{z) and integrate over z to get the 2D 
equation | 22 l EH] 


.d(j)2D{P,t) 
‘ dt 


S + i!f!±rV+ 


y/^dz 


d'k 

-^e~^'^'’Pn{kp,t)h2D 

(27r)^ 


72 


' Mp,tf 


where kp = ^Jkf + ky, and 

n(kp,t)= J dpe*'"'’-^|(/>2D(p,t)7 n{kz) = j dze*''"^|^(z)p = e“77/4. 
1 r°° 

h2n[i)=^J_J 


?>kl 

k2 


^(7)p= d-, [2 - 37^Cexp(g7{l - erf(C)}], f= 


\/^dz 


72 ■ 


(46) 

(47) 

(48) 
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To derive Eq. (46 1 , the dipolar term in Eq. (|^ is first written in momentum space using Eq. 
integration over kz is performed in the dipolar term. 


and the 


The 2D GP equation (46 1 can be solved analytically using the Lagrangian variational formalism with 


the following Gaussian ansatz for the wave function 


7r-i/2 

</>2 d(p,0 = — Try exp 
Wp{t} 


2wj{t) 


+ ia{t)p^ 


(49) 


where Wp(t) is the width and a(t) is the chirp. The Lagrangian variational formalism leads to the following 
equation for the width Wp | 22 ) : 


Wp{t) + j'^Wpit) 


1 A^at [2a - Qddg(K)] _ _ Wp{t) 

Wp(t) ^ dz 


(50) 


The time-independent width of a stationary state can be obtained from Eq. (50) by setting Wp(t) = 0. The 
variational chemical potential for the stationary state is given by [ 22 ] 


1 2jVat[a- add/(K)] 

y/^dzw'^p 2 


(51) 


The energy per particle is given by 

E _ iVat[a - addfjk)] 

Nat “^Wp \/^dzWp 2 


(52) 


2.6.2. X — z plane 

For a disk-shaped dipolar BEG with a strong axial trap along y direction (z/ > A, 7 ), we assume that the 
dynamics of the BEG in the y direction is confined in the ground state 


(t>{y) = exp(-^/V2dy)/(7^c^y)^/^ dy = 07^, 
and we have for the wave function 

(t>{T) = (t>{y) X (t> 2 D{P,t) = 


1 


{ndlYN 


exp 


y 

'2dl, 


4>2d{P, t), 


(53) 


(54) 


where now p = (x,z), and ^ 2 d(p, t) is the circularly-asymmetric effective 2D wave function for the 2D 
dynamics and dy is the harmonic oscillator length along y direction. To derive the effective 2D equation for 
the disk-shaped dipolar BEG, we use ansatz (54) in Eq. ([^, multiply by the ground-state wave function 
4>[y) and integrate over y to get the 2D equation 


.d(j)2D{P,t) 
' dt 


p 


^T:aNat\4‘2D\ 


t/^dy 


d-TTOdd^at 


dk 




{2ny 


'(kp,t)j 


2D 


72 


4^2D{p.t), 


(55) 


where kp = y/k"^ + k'^, 

and 

1 f°° 

Nd{ 0 = ^J dky 

■3kl ■ 

k 2 


\n{ky)? = [- 1 + 377 r^exp(g 7 {l-erf(^)}], ^ 

(56) 

To derive Eq. (55), the dipolar term in Eq. ([^ is first written in momentum space using Eq. ([^ and the 
integration over ky is performed in the dipolar term. 
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3. Details about the programs 


3.1. Description of the programs 

In this subsection we describe the numerical codes for solving the dipolar GP equations (29) and (41) in 
ID, Eq. (46) and (55) in 2D, and Eq. ^ in 3D using real- and imaginary-time propagations. The real-time 
propagation yields the time-dependent dynamical results and the imaginary-time propagation yields the 
time-independent stationary solution for the lowest-energy state for a specific symmetry. We use the split- 
step Crank-Nicolson method for the solution of the equations described in Ref. [3]. The present programs 
have the same structure as in Ref. [3] with added subroutines to calculate the dipolar integrals. In the 
absence of dipolar interaction the present programs will be identical with the previously published ones [3]. 

A general instruction to use these programs in the nondipolar case can be found in Ref. [5] and we refer the 
interested reader to this article for the same. 

The present Fortran programs named (‘imagldX.fQO’, ‘imagldZ.f90’), (‘imag2dXY.f90’, ‘imag2dXZ.f90’), 
‘imag3d.f90’, (‘realldX.f90’, ‘realldZ.f90’), (‘real2dXY.f90’, ‘real2dXZ.f90’), ‘real3d.f90’, deal with imaginary- 
and real-time propagations in ID, 2D, and 3D and are to be contrasted with previously published pro¬ 
grams [3] ‘imagtimeld.F’, ‘imagtime2d.f90’, ‘imagtime3d.f90’, ‘realtimeld.F’, ‘realtime2d.f90’, and ‘real- 
time3d.f90’, for the nondipolar case. The input parameters in Fortran programs are introduced in the 
beginning of each program. The corresponding C codes are called (imagldX.c, imagldX.h, imagldZ.c, 
imagldZ.h,), (imag2dXY.c, imag2dXY.h, imag2dXZ.c, imag2dXZ.h,), (imag3d.c, imag3d.h), (realldX.c, 
realldX.h, realldZ.c, realldZ.h,), (real2dXY.c, real2dXY.h, real2dXZ.c, real2dXZ.h,), (realSd.c, realSd.h), 
with respective input files (‘imagldX-input’, ‘imagldZ-input’), (‘imag2dXY-input’, ‘imag2dXZ-input’), ‘imag3d- 
input’, (‘realldX-input’,‘realldZ-input’), (‘real2dXY-input’, ‘realldXZ-input’), ‘real3d-input’, which per¬ 
form identical executions as in the Fortran programs. 

We present in the following a description of input parameters. The parameters NX, NY, and NZ in 3D 
(NX and NY in 2DXY, NX and NZ in 2DXZ), and N in ID stand for total number of space points in x, y 
and 2 directions, where the respective space steps DX, DY, and DZ can be made equal or different; DT is 
the time step. The parameters NSTP, NPAS, and NRUN denote number of time iterations. The parameters 
GAMMA ( 7 ), NU (ly), and LAMBDA (A) denote the anisotropy of the trap. The number of atoms is denoted 
NATOMS (iVat), the scattering length is denoted AS (a) and dipolar length ADD (odd)- The parameters 
GO (47rA"ata) and GDDO (3addAat) are the contact and dipolar nonlinearities. The parameter OPTION = 2 
(default) defines the equations of the present paper with a factor of half before the kinetic energy and trap; 
OPTION = 1 defines a different set of GP equations without these factors, viz Ref. [3]. The parameter AHO 
is the unit of length and Bohr_aO is the Bohr radius. In ID the parameter DRHO is the radial harmonic 
oscillator dp and in 2D the parameter D_Z or D_Y is the axial harmonic oscillator length dz or dy. The 
parameter CUTOFF is the cut-off R of Eq. (|^ in the 3D programs. The parameters GPAR and GDPAR are 
constants which multiply the nonlinearities GO and GDDO in realtime routines before NRUN time iterations 
to study the dynamics. 

The programs, as supplied, solve the GP equations for specific values of dipolar and contact nonlinearities 
and write the wave function, chemical potential, energy, and root-mean-square (rms) size(s) etc. For solving 
a stationary problem, the imaginary-time programs are far more accurate and should be used. The real-time 
programs should be used for studying non-equilibrium problems reading an initial wave function calculated 
by the imaginary-time program with identical set of parameters (set NSTP = 0, for this purpose, in the 
real-time programs). The real-time programs can also calculate stationary solutions in NSTP time steps (set 
NSTP 7 ^ 0 in real-time programs), however, with less accuracy compared to the imaginary-time programs. 
The larger the value of NSTP in real-time programs, more accurate will be the result [3]. The nonzero 
integer parameter NSTP refers to the number of time iterations during which the nonlinear terms are 
slowly introduced during the time propagation for calculating the wave function. After introducing the 
nonlinearities in NSTP iterations the imaginary-time programs calculate the final result in NPAS plus NRUN 
time steps and write some of the results after NPAS steps to check convergence. The real-time programs 
run the dynamics during NPAS steps with unchanged initial parameters so as to check the stability and 
accuracy of the results. Some of the nonlinearities are then slightly modified after NPAS iterations and the 
small oscillation of the system is studied during NRUN iterations. 
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Each program is preset at fixed values of contact and dipolar nonlinearities as calculated from input 
scattering length(s), dipolar strength(s), and number of atom(s), correlated DX-DT values and NSTP, 
NPAS, and NRUN, etc. A study of the correlated DX and DT values in the nondipolar case can be found 
in Ref. [5]. Smaller the steps DX, DY, DZ and DT, more accurate will be the result, provided we integrate 
over a reasonably large space region by increasing NX, NY, and NZ, etc. Each supplied program produces 
result up to a desired precision consistent with the parameters employed — GO, GDDO, DX, DY, DZ, DT, 
NX, NY, NZ, NSTP, NPAS, and NRUN etc. 

3.2. Description of Output files 

Programs ‘imagnd*’ (n = 1, 2, 3, G and Fortran): They write final density in files ‘imagnd-den.txt’ after 
NRUN iterations. In addition, in 2D and 3D, integrated ID densities ‘imagud*-denldjx.txt’, ‘imagnd*- 
denld_y.txtViniagnd*-denld_z.txt’, along x, y, and z, etc, are given. These densities are obtained by inte¬ 
grating the densities over eliminated space variables. In addition, in 3D integrated 2D densities ‘imagSd- 
den2djcy.txt’, ‘imag3d-den2d_yz.txt’,‘imag3d-den2d_zx.txt’, in xy, yz, and zx planes can be written (com¬ 
mented out by default). The files ‘imagnd*-out.txt’ provide different initial input data, as well as chemical 
potential, energy, size, etc at different stages (initial, after NSTP, NPAS, and after NRUN iterations), from 
which a convergence of the result can be inferred. The files ‘imagnd*-rms.txt’ provide the different rms sizes 
at different stages (initial, after NSTP, NPAS, and after NRUN iterations). 

Programs ‘realnd*’ (n = 1, 2,3, C and Fortran): The same output files as in the imaginary-time programs 
are available in the real-time programs. The real-time densities are reported after NPAS iterations. In 
addition in the ‘realud*-dyna.txt’ file the temporal evolution of the widths are given during NPAS and 
NRUN iterations. Before NRUN iterations the nonlinearities GO and GDDO are multiplied by parameters 
GPAR and GDPAR to start an oscillation dynamics. 

3.3. Running the programs 

In addition to installing the respective Fortran and C compilers one needs also to install the EFT routine 
FFTW in the computer. To run the Fortran programs the supplied routine fftw3.f03 should be included in 
compilation. The commands for running the Fortran programs using INTEL, GFortran, and Oracle Sun 
compilers are given inside the Fortran programs. The programs are submitted in directories with option 
to compile using the command ‘make’. There are two files with general information about the programs 
and FT for user named ’readme.txt’ and ‘readme-fftw.txt’. The Fortran and G programs are in directories 
../Lprogram and ../c.program. Inside these directories there are subdirectories such as ../input, ../output, 
../src. The subdirectory ../output contains output files the programs generate, ../input contains input 
files for C programs, and ../src contains the different programs. The command ‘make’ in the directory 
../Lprogram or ../c-program compiles all the programs and generates the corresponding executable files to 
run. The command ‘make’ for INTEL, GFortran and OracleSun Fortran are given. 


4. Numerical Results 


In this section we present results for energy, chemical potential and root-mean-square (rms) sizes for 
different stationary BEGs in ID, 2D, and 3D, and compare with those obtained by using Gaussian variational 
and TF approximations, wherever possible. We also compare with available results by other authors. For a 
fixed space and time step, sufficient number of space discretizing points and time iterations are to be allowed 
to get convergence. 

First we present in table numerical results for the energy E, chemical potential /i, and rms size (z) 
calculated using the imaginary-time program for the ID dipolar GP equation (29) for ®^Gr atoms with 
0 = 6 nm (« llSoo with oq the Bohr radius), and Odd = 16ao for A = 1, dp = 1, Z = 1 /im and for different 
number of atoms Nat and different space and time steps dz and dt. The Gaussian variational approximations 


obtained from Eqs. (36), (37) and (38) are also given for comparison. The variational results provide better 


approximation to the numerical solution for a smaller number of atoms. 
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Table 1: The energy per particle E/N^t, chemical potential /i., and rms size (z) of the ID GP equation ( |29[> f or X = 1, dp = 1 
fim for the ^^Cr BEG with a = 6 nm, Udd = 16ao and different number of atoms A^at- In Eqs. and ( |29[ l the lengths are 
expressed in oscillator unit: I = 1 fim. Numerical results are calculated for parameters (A) dz = 0.05, dt = 0.0005, A = 2048 
(B) dz = 0.1, dt = 0.001, N = 1024 and compared with variational results obtained from Eqs. ( |36[ l, l |37[ | and | |38[ |. 


iVat 

(^) 

var 

{z) 

(B) 

{z) 

(A) 

A/iVat 

var 

E/N,,t 

(B) 

E/N^t 

(A) 

var 

d 

(B) 

(A) 

100 

0.7939 

0.7937 

0.7937 

0.7239 

0.7222 

0.7222 

0.9344 

0.9297 

0.9297 

500 

1.0425 

1.0381 

1.0381 

1.4371 

1.4166 

1.4166 

2.2157 

2.1691 

2.1691 

1000 

1.2477 

1.2375 

1.2375 

2.1376 

2.0920 

2.0920 

3.4165 

3.3234 

3.3234 

5000 

2.0249 

1.9939 

1.9939 

5.8739 

5.6910 

5.6910 

9.6671 

9.3488 

9.3488 

10000 

2.5233 

2.4815 

2.4815 

9.2129 

8.913 

8.913 

15.223 

14.715 

14.715 

50000 

4.2451 

4.1719 

4.1719 

26.505 

25.622 

25.622 

43.993 

42.527 

42.527 


Table 2: The energy per particle E/Nat, chemical potential fi, and rms size (p) of the 2D GP equation ( |46| l for '7 = u = 1, dz = 1 
pm for the ^^Cr BEG with a = 6 nm, a^d = 16ao and different number of atoms N^f In Eqs. and | |46| the lengths are 
expressed in oscillator unit: I = 1 /im. Numerical results are calculated for space and time steps (A) dx = dy = 0.1, dt = 
0.0005, NX = NY = M = 768, (B) dx = dy = 0.2, dt = 0.002, A7 = 384, and compared with variational results obtained from 
Eqs. l |50| l, ( |51| l and ( |52[ |. 


A^at 

(P) 

var 

(P) 

(B) 

(P) 

(A) 

E/N^t 

var 

A/iVat 

(B) 

A/Vat 

(A) 

var 

P 

(B) 

P 

(A) 

100 

1.0985 

1.097 

1.097 

1.2182 

1.2156 

1.2157 

1.4187 

1.4120 

1.4119 

500 

1.3514 

1.342 

1.342 

1.8653 

1.8383 

1.8383 

2.5437 

2.4840 

2.4840 

1000 

1.5482 

1.530 

1.531 

2.4571 

2.3988 

2.3988 

3.5070 

3.3901 

3.3901 

5000 

2.2549 

2.208 

2.208 

5.2206 

4.9989 

4.9989 

7.8005 

7.4249 

7.4249 

10000 

2.6824 

2.619 

2.619 

7.3787 

7.029 

7.029 

11.090 

10.522 

10.522 

50000 

4.0420 

3.934 

3.934 

16.680 

15.793 

15.793 

25.161 

23.789 

23.789 


In tab le we present results for the energy E, chemical potential /i, and rms size (p) of the 2D GP 
equation ( |46[ ) for 7 = t/ = l,dz = 1,1 = 1 /im. The numerical results are calculated using different space 
and time steps dx, dy and dt and different number TVat of ®^Cr atoms with Odd = ISoq and a = 6 nm. 
Axially-symmetric Gaussian variational approximations obtained from Eqs. (50), (51) and (52) are also 
presented for comparison. 

Now we present results of the solution of the 3D GP equation (|^ with some axially-symmetric traps. In 
this case we take advantage of the cut-off introduced in Eq. to improve the accuracy of the numerical 
calculation. The cut-off parameter R was taken larger than the condensate size and smaller than the 
discretization box. First we consider the model 3D GP equation with a = 0 and different ^dd = = 

1,2, 3,4 in an axially-symmetric trap with A = 1/2 and z/ = 7 = 1. The numerical results for different 
number of space and time steps together with Gaussian variational results obtained from Eqs. (18) and 
(19) are shown in table These results for energy E and chemical potential p are compared with those 
calculated by Asad-uz-Zaman et al. [I6l[29]. The present calculation is performed in the Gartesian x,y,z 
coordinates and the dipolar term is evaluated by FT to momentum space. Asad-uz-Zaman et al. take 
advantage of the axial symmetry and perform the calculation in the axial p,z {p = x,y) variables and 
evaluate the dipolar term by a combined Hankel-Fourier transformation to momentum space for p and z, 
respectively. The calculations of Asad-uz-Zaman et al. for stationary states involving two variables {p and z) 
thus could be more economic and accurate than the present calculation involving three Gartesian variables 
for the axially-symmetric configuration considered in table However, the present method, unlike that of 
Ref. [12], is readily applicable to the fully asymmetric configurations. Moreover, the present calculation for 
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Table 3: Energy per particle E/Ns^t and chemical potential p from a solution of Eq. § for 7 = z/ = 1, = 0.25, a = 0 and 

different nonlinearity = 3add^at- The present numerical results are compared with Gaussian variational results obtained 
from Eqs. ( |18[ | and l |19[ | as well as numerical results of Asad-uz-Zaman et al. [la I29| . Numerical results are calculated 
for the following space and time steps and the following space discretizing points in the Crank-Nicolson discretization: (A) 
dx = dy = dz = 0.05, dt = 0.0004, {NX = NY = NZ = M = 384); (B) 0.1, dt = 0.002, {N = 128, = 6); and (C) 

0.2, dt = 0.007, {N = 64, R = 6). 


9dd 

E/N,t 

var 

E/N,,t 

(C) 

E/N,t 

(B) 

E/N,t 

(A) 

E/N,,t 

m 

R 

var 

(C) 

(B) 

(A) 

m 

0 

1.2500 

1.2498 

1.2500 

1.2500 

1.2500 

1.2500 

1.2498 

1.2500 

1.2500 

1.2500 

1 

1.2230 

1.2220 

1.2222 

1.2222 

1.2222 

1.1934 

1.1910 

1.1912 

1.1911 

1.1911 

2 

1.1907 

1.1872 

1.1875 

1.1874 

1.1874 

1.1203 

1.1100 

1.1100 

1.1100 

1.1100 

3 

1.1521 

1.143 

1.1439 

1.1438 

1.1437 

1.0253 

0.995 

0.996 

0.996 

0.9955 

4 

1.1051 

1.085 

1.0857 

1.0857 

1.0856 

0.8950 

0.805 

0.803 

0.806 

0.8062 


dynamics (non-stationary states) in 3D are more realistic than the calculations of Asad-uz-Zaman et al., 
where one degree of freedom is frozen. For edxample, a votex could be unstable |30j in a full 3D calculation, 
whereas a 2D calculation could make the same vortex stable. 


Table 4: Energy per particle E/N^t, and chemical potential fi from a solution of Eq. § for 7 = zz = 1, = 0.25, Ana = 

0.20716, 47radd = 0.033146 and different number Nat of atoms. These nonlinearity parameters taken from Ref. [9] correspond 
to a ^^Cr dipolar BEG with a ps IOOuq and a^d ~ 16ao and oscillator length I ps 0.321 pm. Variational and TF results as 
well as numerical results of Bao et al. [9] are also shown. Numerical results are calculated using the following space and time 
steps and the following space discretizing points in the Crank-Nicolson discretization: (A) dx = dy = dz = 0.15, dt = 0.002; 
(B) dx = dy = dz = 0.3, di = 0.005. In (A) we take NX = NY = NZ = J\f = 128, R = 9 for Nat = 100,500,1000 
and J\f = 192, i? = 14, for Nat = 5000,10000,50000 ; and in (B) we take Af = 64, i? = 9, for Nat = 100,500,1000 and 
Af = 96,R= 14, for Vat = 5000,10000, 50000. 


N^t 

E/N,t 

var 

E/N,t 

TF 

E/N,t 

(B) 

E/N,,t 

(A) 

E/N,t 

m 

R 

var 

IJ- 

TF 

(B) 

(A) 

m 

100 

1.579 

0.945 

1.567 

1.567 

1.567 

1.840 

1.322 

1.813 

1.813 

1.813 

500 

2.287 

1.798 

2.224 

2.224 

2.225 

2.951 

2.518 

2.835 

2.835 

2.837 

1000 

2.836 

2.373 

2.728 

2.728 

2.728 

3.767 

3.322 

3.583 

3.582 

3.583 

5000 

5.036 

4.517 

4.744 

4.744 

4.745 

6.935 

6.324 

6.485 

6.486 

6.488 

10000 

6.563 

5.960 

6.146 

6.146 

6.147 

9.100 

8.344 

8.475 

8.475 

8.479 

50000 

12.34 

11.35 

11.46 

11.46 

11.47 

17.23 

15.89 

15.96 

15.97 

15.98 


Table 5: The rms sizes (x) and (z) for the same systems illustrated in 


table using the same cut-off parameter R. 


N {z) {z) {z) {z) (z) {x) (x) (x) (x) (x) 



TF 

var 

(B) 

(A) 

m 

TF 

var 

(B) 

(A) 

m 

100 

1.285 

1.316 

1.305 

1.303 

1.299 

0.600 

0.799 

0.794 

0.795 

0.796 

500 

1.773 

1.797 

1.752 

1.752 

1.745 

0.828 

0.952 

0.938 

0.939 

0.940 

1000 

2.037 

2.079 

2.014 

2.014 

2.009 

0.951 

1.054 

1.035 

1.035 

1.035 

5000 

2.810 

2.904 

2.795 

2.795 

2.790 

1.313 

1.392 

1.353 

1.353 

1.354 

10000 

3.228 

3.345 

3.217 

3.216 

3.212 

1.508 

1.586 

1.537 

1.537 

1.538 

50000 

4.454 

4.629 

4.450 

4.450 

4.441 

2.080 

2.171 

2.093 

2.093 

2.095 
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Table 6; Energy per particle E/Naf, chemical potential p, and rms sizes from a solution of Eq. ^ for ^^Cr atoms with 
7 = 1,1^^ = 1/2, = 1/4, a = llOao, CLdd = 16ao, and harmonic oscillator length I = 1 pm for different A^at- Numerical 

results are calculated using the following space and time steps and the following space discretizing points in the Crank-Nicolson 
discretization: (A) dx = dy = dz = 0.1, dt = 0.001; and (B) 0.2, dt = 0.003. In (A) we take NX = NY = NZ = J\f = 128, R = 6 
for A^at = 100,500,1000 and J\f = 256, R = 10.5, for N^t = 5000,10000,50000 ; and in (B) we take J\f = 64, i? = 6, for 
Aat = 100, 500,1000 and Af = 128, R = 12, for Aat = 5000,10000, 50000. 


N 

E/N,t 

(B) 

E/N,t 

(A) 

d 

(B) 

d- 

(A) 

(x) 

(B) 

(y) 

(B) 

(z) 

(B) 

(x) 

(A) 

(y) 

(A) 

(^) 

(A) 

100 

1.219 

1.219 

1.321 

1.321 

0.742 

0.901 

1.120 

0.742 

0.901 

1.119 

500 

1.525 

1.525 

1.830 

1.830 

0.818 

1.032 

1.379 

0.818 

1.032 

1.379 

1000 

1.784 

1.784 

2.232 

2.232 

0.874 

1.128 

1.559 

0.874 

1.129 

1.558 

5000 

2.885 

2.885 

3.857 

3.858 

1.079 

1.463 

2.132 

1.079 

1.463 

2.132 

10000 

3.673 

3.673 

4.992 

4.992 

1.206 

1.660 

2.450 

1.206 

1.660 

2.449 

50000 

6.713 

6.713 

9.306 

9.306 

1.609 

2.260 

3.383 

1.609 

2.260 

3.383 


Next we consider the solution of the 3D GP equation ([^ for a model condensate of ®^Cr atoms in a cigar¬ 
shaped axially-symmetric trap with 7 = z/ = l,A = l/2, first considered by Bao et al. [3. The nonlinearities 
considered there (47ra = 0.20716, dTrodd = 0.033146) correspond to the following approximate values of a, Odd 
and 1: a « IOOuq, add ~ 16ao, and I = 0.321 /rm. We present results for energy E and chemical potential ^ 
in table 1^ and rms sizes {z) and {x) in table We also present variational and Thomas-Fermi (TF) results 
in this case together with results of numerical calculation of Bao et al. [3]. The TF energy and chemical 
potential in tableare calculated using Eqs. ( [^ and (26), respectively. The TF sizes {x) and {z) in table 
are obtained from Eqs. (24) and (25) using the TE density (20). For small nonlinearities or small number of 


atoms, the Gaussian variational results obtained from Eqs. (13), (14), (18), and (19) are in good agreement 


with the numerical calculations as the wave function for small nonlinearities has a quasi-Gaussian shape. 
However, for large nonlinearities or large number of atoms, the wave function has an approximate TF shape 
(^, and the TF results provide better approximation to the numerical results, as can be seen from tables 

^andlll 

After the consideration of 3D axially-symmetric trap now we consider a fully anisotropic trap in 3D. In 
table 1^ we present the results for energy E/N^t, chemical potential /r and rms sizes (x), (y), (z) of a ®^Cr 
BEG in a fully anisotropic trap with 7 = 1,1/= l/-\/2, A = 1/2 for different number of atoms. In this case we 
take a = llOoo, Odd = 16ao and I = 1 /rm. The convergence of the calculation is studied by taking reduced 
space and time steps dx and dt and different number of space discretization points. Sufficient number of 
time iterations are to be allowed in each case to obtain convergence. In 3D the estimated numerical error 
in the calculation is less than 0.05%. The error is associated with the intrinsic accuracy of the FFT routine 
for long-range dipolar interaction. 


The ID and 2D GP equations (29) and (46) are valid for cigar- and disk-shaped BEGs, respectively. In 


case of cigar shape the ID GP equation yields results for axial density and in this case it is appropriate 
to compare this density with the reduced axial density obtained by integrating the 3D density over radial 

1 (a) we compare two axial densities obtained 


coordinates: n{z) = \^p{z)'\^ = J \(l){x, y, z)f dx dy. In Fig. 


from the ID and 3D GP equations. We also show the den^ies calculated from the Gaussian variational 
approximation in both cases. In the cigar case the trap parameters are i/ = 7 = 1,A = 1/4. Similarly, for 
the disk shape it is interesting to compare the density along the radial direction in the plane of the disk as 
obtained from the 3D equation ([^ and the 2D equation (461. In this case it is appropriate to calculate the ID 
radial density along, say, x direction by integrating 2D and 3D densities as follows: nioix) = f dy\4>2Dix, j/)p 
and nioix) = f dydz\(l) 3 D{x,y, z)\‘^ . In Fig. 0 (b) we compare two radial densities obtained from the 2D 
and 3D GP equations. We also show the densities calculated from the Gaussian variational approximation 
in both cases. For this illustration, we consider the trap parameters z/ = 7 = l,A = 4. In both Figs. (a) 


15 



























Figure 1: (a) Numerical (num) and variational (var) results for the one-dimensional axial density 711 / 3 ( 2 :) = | 0 i/ 3 ( 2 )p along 2 
axis for 1 / = 7 = 1, A = 0.25 of a cigar-sh^ed BEC of N^t = 1000 atoms obtained using the ID Eq. ( |29[ l and that obtained 
after integrating the 3D density from Eq. lliMI over x and y. 711 / 3 ( 2 ) = f \(l>{r)\^dxdy. (b) Numerical (num) and variational (var) 
results for the ID radial density ni/ 3 (x) = J \(l>{r)\^dydz along x axis for M = 7 = l,A = 4ofa disk-shaped BEC of iVat = 1000 
atoms obtained after integrating the 3D density from Eq. over y and 2 and after integrating the 2D density from Eq. l |46[ l 
over y as follows: 711 / 3 ( 2 ;) = f dy\<p 2 D{^ty)\‘^ a-nd = J d7/(i2|03/3(2;, j/, 2 )^. In all cases a = 6 nm and = 16ao. 




Figure 2: (a) The Nat — o. stability phase plot for a ^®^Dy BEC with a,/;! = 130ao in a disk-shaped trap with i/ = 7 = l, A = 5 
and 7 and harmonic oscillator length I = 1 /im. The 3D isodensity contour plot of density of a disk-shaped ^®'*Dy BEC with 
Ujiu = 130ao for 71 = 7 = 1 , A = 5, i = l /im, Nat = 3000 and a = 40ao for densities \(l>{x, y, 2 )^ = (b) 0.001 and (c) 0.027 on 
the contour. 


and (b), the densities obtained from the solution of the 3D GP equation are in satisfactory agreement with 
those obtained from a solution of the reduced ID and 2D equations. In Fig. 1, the numerical and variational 
densities are pretty close to each other, so are the results obtained from the 3D equation ([^, on the one 
hand, and the ones obtained from the ID and 2D equations (29) and (46), on the other. 

A dipolar BEC is stable for the number of atoms TVat below a critical value [3T]. Independent of trap 
parameters, such a BEC collapses as iVat crosses the critical value. This can be studied by solving the 3D 
GP equation using imaginary-time propagation with a nonzero value of NSTP while the nonlinearities are 
slowly increased. In Pig. (a) we present the N^t — a stability phase plot for a ^^‘^Dy BEC with add = 130ao 
in the disk-shaped trap with v = ^ = l, X = 5 and 7. The oscillator length is taken to be ( = I ^m. The 
shaded area in these plots shows a metastable region where biconcave structure in 3D density appears. The 
metastable region corresponds to a local minimum in energy in contrast to a global minimum for a stable 
state. It has been established that this metastability is a manifestation of roton instability encountered by 
the system in the shaded region |31j . The biconcave structure in 3D density in a disk-shaped dipolar BEC 
is a direct consequence of dipolar interaction: the dipolar repulsion in the plane of the disk removes the 
atoms from the center to the peripheral region thus creating a biconcave shape in density. In Figs. Kb) 
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Figure 3: (a) Numerical results for the ID radial density ni£)(x) = f \4){r)\’^dydz along x axis and ni£)( 2 :) = f \(l){r)\^dxdy 

along 2 axis for A = 7 = l,z/ = 4ofa disk-shaped BEC of Nat = 2000 atoms obtained after integrating the 3D density 

from Eq. and the 2D density from Eq. ( |55[ ) over the eliminated variables, (b) Numerical results for the ID axial den sity 
^id(^) = along X axis for u = \ = 16, 7 = 1 of a cigar-shaped BEC of Nat = 2000 ^^Cr atoms obtained using the ID Eq. |4y 
and that obtained after integrating the 3D density from Eq. over z and y: ni£){x) = f \(l){r)\‘^dzdy. In all cases a = 120ao 
and (a) add = 132.7ao, (b) add = 16ao. 


and (c) we plot the 3D isodensity contour of the condensate for A — 5 with parameters in the shaded region 
corresponding to metastability. In Fig. (b) the density on the contour is 0.001 whereas in Fig. Kc), 
it is 0.027. Only for a larger density on the contour the biconcave shape is visible. The biconcave shape 
predominates near the central region of the metastable dipolar BEC. 

In figure we critically tested the reduced ID and 2D equations (291 and (46) along the z axis and 
in the x — y plane, respectively, by comparing the different ID densities from these equations with those 
obtained from a solution of the 3D equation as well as with the variational densities. Now we perform 
a similar test with the reduced ID and 2D equations (41) and (55) along the x axis and in the x — z plane, 
respectively. We consider a BEC of 2000 atoms in a disk-shaped trap in the x — z plane with A = 7 = 1 and 
V = A. Because of the strong trap in the y direction, the resultant BEC is of quasi-2D shape in the x — z 
plane without circular symmetry in that plane because of the anisotropic dipolar interaction. The integrated 
linear density along the x and z axes as calculated from the 2D GP equation (55) and the 3D GP equation 
(§ are illustrated in figure]^ (a). Next we consider the BEC of 2000 atoms in a cigar-shaped trap along the 
X axis with v = X = IQ and 7 = 1 . The integrated linear density along the x axis in this case calculated 
from the 3D equation ([^ is compared with the same as calculate using the reduced ID equation (41) in 
figure 1^ (b). In both cases the densities calculated from the 3D GP equation are in reasonable agreement 
with those calculated using the reduced equations (55) and ( |4T| ). Another interesting feature emerges from 
figures and 1^ the reduced 2D GP equations (46) and ( [^ with appropriate disk-shaped traps yield results 
for densities in better agreement with the 3D GP equation ^ as compared to the ID GP equations ( [^ 
and ( [ 4 I] ) with appropriate cigar-shaped traps. This feature, also observed in non-dipolar BECs [T7], is 
expected as the derivation of the reduced ID equations involving two spatial integrations represent more 
drastic approximation compared to the same of the reduced 2D equations involving one spatial integration. 

Now we report the dynamics of the dipolar BEC by real-time propagation using the stationary state 
calculated by imaginary-time propagation. In Fig . [4] (a) we show the oscillation of the rms sizes (z) and 
(p) from the reduced ID and 2D GP equations ( |29[ ) and ( [4^ , respectively. In Fig. (a) we consider 
A^at = 10000, Odd = ISoo (appropriate for ®^Cr), a = 6 nm (« I13ao) and oscillator length I = I ^m. In 
ID, we take dx = 0.025, dt = 0.0001, A = 1, dp = 1, number of space points N = 2048, and in 2D, we take 
dx = dy = 0.2, dt = 0.001,7 = v = l,dz = \,NX = NY = 512.in real-time simulation the oscillation is 
started by multiplying the nonliniarities with the factor 1.05. To impliment this, in real-time routine we take 
GPAR = GDPAR = 1.1 and also take NSTP = 0 to read the initial wave function. In ID and 2D we also 
present results of the Gaussian variational approximations after a numerical solution of Eqs. (36) and (50), 
respectively. The frequency of the resultant oscillations agree well with the numerical ID and 2D calculations. 
However, slight adjustment of the initial conditions, or initial values of width and its derivative, were 
necessary to get an agreement of the amplitude of oscillation obtained from variational approximation and 
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Figure 4: (a) Numerical (n) and variational (u) results for oscillation of rms sizes {z) and (p) from the real-time simulation 
using Eq. | |29[ | in ID and Eq. | |46[ l in 2D, respectively, for A^at = 10000, a = 6 nm, a^d = lOo-O) ^ = 1 while a and were 
both multiplied by 1.05 after NPAS iterations at £ = 10. The wave function was first calculated by imaginary-time routine with 
parameters dx = 0.025, dt = 0.0001, A = l^dp = 1, NPAS= 10^, A" = 2048 in ID, and dx = dy = 0.2, dt = 0.001,7 = 1,^2 = 1 
NPAS = 10“^, NX = NY = 512 i n 2 D. The results of the variational approximations in ID and 2D as obtained from a 
numerical solution of Eqs. ( |36[ ) and ( |50[ l are also shown, (b) Numerical results for oscillation of rms sizes (a;), {y) and {z) from 
the real-time simulation in 3D using Eq. §, for Aat = 1000, a = llOap, ^dd = lOctO) ^ = \^u = l/\/2;A = 1/2, 

NX = NY = NZ = 128, da: = dy = dz = 0.2, and dt = 0.002 while a and add were both multiplied by 1.1 after NPAS 
iterations. In all cases the real-time calculation was performed with NSTP = 0 reading the 3D density from the numerical 
solution of the imaginary-time program using the same parameters. 


numerical simulation. The initial values of width and its derivative are necessary to solve the variational Eqs. 
(36l and (50). In Fig. |^(b) we illustrate the oscillation of the rms sizes (x), (y), and (z) in 3D using Eq. (§, 
where we perform real-time simulation using the bound state obtained by imaginary-time simulation as the 
initial state. The parameters used are N^t = 1000,a = 110ao,add = 16ao)7 = l,v = 11^/2,\ = 1/2,1 = 1 
/im, NX = NY = NZ = 128, dx = dy = dz = 0.2, dt = 0.002 in both real- and imaginary-time simulations. 
In addition, in real-time simulation the oscillation is started by multiplying the nonliniarities with the factor 
1.1. To impliment this, in real-time routine we take GPAR = GDPAR = 1.1 and also take NSTP = 0 to 
read the initial wave function. 


5. Summary 

We have presented useful numerical programs in Fortran and C for solving the dipolar GP equation 
including the contact interaction in ID, 2D, 3D. Two sets of programs are provided. The imaginary- 
time programs are appropriate for solving the stationary problems, while the real-time codes can be used 
for studying non-stationary dynamics. The programs are developed in Cartesian coordinates. We have 
compared the results of numerical calculations for statics and dynamics of dipolar BECs with those of 
Gaussian variational approximation, Thomas-Fermi approximation, and numerical calculations of other 
authors, where possible. 
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